Modul 12 - Contoh Analisis Big Data

Pengantar

Data Time Series adalah data yang diperoleh dari pengamatan satu objek dari beberapa periode waktu. Misalnya data jumlah mahasiswa dari tahun ke tahun, data nilai tukar dolar terhadap rupiah, data kasus covid dari bulan januari 2020 sampai periode saat ini, dan masih banyak contoh yang lainnya. Analisis Time Series adalah suatu bentuk peramalan terhadap nilai-nilai dimasa yang akan datang yang didasarkan pada nilai-nilai pada masa lampau. Model ini biasanya digunakan untuk melakukan prediksi/peramalan. Untuk melakukan analisis time series tersebut ada banyak teknik/metode/algoritma/model yang dapat digunakan, berikut beberapa diantaranya:

1. ARIMA, digunakan untuk tipe data stationary / random
2. ETS (error, trend, seasonal), sangat mirip pendekatannya dengan prophet facebook, yaitu dengan komponen trend dan seasonality
3. Bayesian Structural Time Series (BSTS), baik digunakan ketika ada eksternal regressor (atau pengaruh eksternal dari series lain selain komponen time series itu sendiri)
4. TBATS, model yang bisa digunakan ketika memodelkan “multiple seasonality” yaitu penggunaan 2 macam efek seasonality (misal weekly dan monthly), karena secara umum time series dianggap hanya memiliki 1 tipe seasonality
5. Prophet adalah model yang dikembangkan oleh Facebook dan fokus pada komponen: trend, seasonality (bisa multiple), dan holiday effect. dll.

Secara umum ada 3 komponen time series yaitu:

• Seosonal : data musiman, data dengan pola yang berulang dengan periode yang tetap (tahun, bulan, minggu, hari, jam, detik)
• Trend : ada 2 jenis trend yaitu cendrung naik atau cendrung turun
• Error : data yang ditangkap oleh Seosonal atau trend yang tidak bisa dijelaskan.

Load Packages

# Library untuk manipulasi data
library (dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

# Library untuk visualisasi data
library (ggplot2)
library (plotly)

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

# Library model time series forecasting
library (prophet)

## Loading required package: Rcpp

## Loading required package: rlang

# Library untuk split data
library (caTools)
# Library manipulasi tanggal
library (lubridate)

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

Data Preparation

Import Data

sales <- read.csv("data/train.csv")
head(sales)

Melihat struktur data

# melihat struktur dataset
glimpse (sales)

## Rows: 730,500
## Columns: 4
## $ date  <chr> "2013-01-01", "2013-01-02", "2013-01-03", "2013-01-04", "2013-01~
## $ store <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1~
## $ item  <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1~
## $ sales <int> 13, 11, 14, 13, 10, 12, 10, 9, 12, 9, 9, 7, 10, 12, 5, 7, 16, 7,~

Exploratory Data Analysis

Menampilkan trankasi berdasarkan store

group_sales <- group_by(sales, store)
count(group_sales)

Bisa dilihat bahwa setiap store memiliki 73050 traksaksi.

Menampilkan jumlah transaksi berdasarkan tanggal

Sekarang kita akan melihat jumlah transaksi berdasarkan tanggal transaksi dan kita simpan dalam objek daily_demand

daily_demand <- sales %>% 
  group_by(date) %>% 
  summarise( 
    demand=sum(sales)
  )  
daily_demand

Menampilkan transaksi pada salah satu store

Kita akan mengambil trankasi pada salah satu store saja, misalnya pada store 3. dan kita simpan dalam objek daily_demandstore03 dan rubah data dateke type tanggal (date)

daily_demandstore03 <- sales %>% 
  filter (store == "3") %>%
  mutate (date = as.Date(date)) %>% #merubah type data date
  group_by(date) %>% 
  summarise( 
    demand=sum(sales)
  ) 
  
daily_demandstore03

Bisa dilihat bahwa data date sudah dalam type tanggal

Visualisasi Data

kita akan memvisualisasi trankasi pada store 03 dari tahun 2013 sampai 2016

plotstore03 <- daily_demandstore03 %>% 
  ggplot(aes(x = date, y = demand)) +
  geom_point(color = "tomato3", group=1) + 
  labs( 
    title = "Daily Sales", 
    subtitle = "Store 03", 
    caption = "1C Company", 
    x = "Date", 
    y = "Total Sales" 
  ) + 
  theme_minimal() 

ggplotly(plotstore03)

Membuat Objek Time Series dengan Model Prophet

Data yang digunakan yaitu data permintaan harian di Store 3 yang sudah disimpan pada objek daily_demandstore03

daily_demandstore03

Menyiapkan data

Untuk menggunakan algoritma/model Prophet pertama-tama yang harus dilakukan adalah menyiapkan data frame dengan format: ds untuk menyimpan tanggal, dan y untuk menyimpan dari nilai yang akan prediksi.

# Menyiapkan data
train_daily_3 <- daily_demandstore03 %>% 
  rename(
    ds = "date",
    y = "demand"
    )
glimpse(train_daily_3)

## Rows: 1,461
## Columns: 2
## $ ds <date> 2013-01-01, 2013-01-02, 2013-01-03, 2013-01-04, 2013-01-05, 2013-0~
## $ y  <int> 1588, 1538, 1635, 1741, 1887, 1956, 1313, 1538, 1633, 1677, 1831, 1~

Fit Model Data Training

Melakukan pemodelan data training dengan menggunakan fungsi fit.prophet() dimana sebelumnya kita harus menentukan seasonality nya. dalam kasus ini kita mengggunakan yearly_seasonality

# fitting model data training
model_ts <- prophet(yearly_seasonality = TRUE, seasonality_prior_scale=0.1) %>% 
  fit.prophet(train_daily_3)

## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.

Berdasarkan hasil fitting model diatas bisa dijelaslkan bahwa model_ts sekarang menyimpan informasi data permintaan harian pada store 03. Sehingga kita bisa gunakan untuk mengekstrak informasi deret waktu tersebut untuk melakukan perkiraan selama periode waktu yang bisa ditentukan. Misalnya kita akan memperkiraan permintaan penjualan untuk 1 tahun kedepan. Utuk melakukan ini kita harus menyiapkan data frame yang terdiri dari rentang waktu/tanggal mendatang yang akan diperkirakan. Untungnya dalam prophet() sudah menyediakan sebuah fungsi make_future_dataframe() memudahkan kita untuk menyiapkan data tersebut.

range(train_daily_3$ds)

## [1] "2013-01-01" "2016-12-31"

Bisa dilihat bahwa jangka waktu prediksi mulai dari 2013 s/d 2016

Menyiapkan tanggal untuk prediksi

Periode prediksi adalah 1 tahun kedepan atau 365 hari

# menyiapkan tanggal utk prediksi
future_ts <- make_future_dataframe(model_ts, periods = 365, freq = "day")
glimpse(future_ts)

## Rows: 1,826
## Columns: 1
## $ ds <dttm> 2013-01-01, 2013-01-02, 2013-01-03, 2013-01-04, 2013-01-05, 2013-0~

Visualisasi hasil forecasting untuk perkiraan data 1 tahun kedepan

# visualisasi hasil forecasting
forecast_ts <- predict(model_ts, future_ts)
plot(model_ts, forecast_ts)

Bisa dilihat bahwa, gambar pada grafik yang berwarna biru menunjukkan hasil forecast yang sudah dilakukan. Titik hitam menunjukkan data dari mulai 2013 sampai 2016. Waktu biru yang tidak ada titik adalah perkiraan permintaan untuk 2017

Visualisasi komponen model

Dengan prophet dimungkinkan kita memvisulisasikan berdasarkan trend, weekly, dan yearly. Default: - daily: sampling hourly, dan frequency 24 - weekly: sampling daily, dan frequency 7 - yearly: sampling daily dan frequency 365

#visualiasi komponen model
prophet_plot_components(model_ts, forecast_ts)

Berdasarkan hasil visualisasi diatas bisa disimpulkan sebagai berikut:

• Penjualan dari tahun ke tahun meningkat, dimana penjualan di tahun 2017 adalah penjualan yang tertinggi
• Penjualan terbanyak pada setiap bulan dan tahun cendrung mengalami peningkatan di hari weekend, yaitu hari Jumat, Sabtu dan Minggu, sementara penjualan di hari Senin adalah penjualan terendah dibandingakan dengn hari yang lain
• Periode penjualan dalam setahun yang paling tinggi adalah pada bulan Juni, Juli dan Agustus, sedangkan tingkat penjualan terendah pada bulan Januari

forecast berdasarkan seasionality bulan

# forecast berdasarkan seasionality month (bulan)
model_ts_monthly <- prophet(changepoint.prior.scale = 0.05, 
                    yearly.seasonality = FALSE) %>% 
  add_seasonality(name = "monthly", period = 30.5, fourier.order = 3) %>% 
  fit.prophet(train_daily_3) 

## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.

future_ts_monthly <- make_future_dataframe(model_ts_monthly, periods = 365) 
forecast_ts_monthly <- predict(model_ts_monthly, future_ts_monthly) 
prophet_plot_components(model_ts_monthly, forecast_ts_monthly)

Visualisasi trend permintaan

Dalam kasu sini kita akan memvisualisasi tren permintaan dari awal tahun 2017 sampai akhir 2017

#visualisasi dari tren awal 2017 sampai akhir 2017
plot(model_ts_monthly, forecast_ts_monthly)

Bisa dilihat bahwa terjadi nya trend peningkatan permintaan yang relatif tinggi dibandingkan tahun sebelumnya

forecast berdasarkan seasionality penambahan holiday (hari libur)

# Menetapkan hari libur mulai dari tahun 2013 sampai 2016, misal libur di hari natal setiap tanggal 25 desember
holiday <- 
  data.frame( 
    holiday = "newyeareve", 
    ds = dmy(c("31-12-2013","31-12-2014", "31-12-2015", "31-12-2016", "31-12-2017")), 
    lower_window = -5, 
    upper_window = 0 
  ) 
holiday 

# visualisasi hasil forecast dengan penambahan efek holiday
model_ts_holiday <- prophet(changepoint.prior.scale = 0.05, 
                    holidays = holiday) %>% 
  add_seasonality(name = "monthly", period = 30.5, fourier.order = 5) %>% 
  fit.prophet(train_daily_3) 

## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.

future_ts_holiday <- make_future_dataframe(model_ts_holiday, periods = 365) 
forecast_ts_holiday <- predict(model_ts_holiday, future_ts_holiday) 
plot(model_ts_holiday, forecast_ts_holiday)

# melihat trend dan holiday effect dari perkiraan
prophet_plot_components(model_ts_holiday, forecast_ts_holiday)

Berdasarkan hasil visualisasi diatas bisa disimpulkan sebagai berikut:

1. Berdasarkan penjualan mingguan didapatkan bahwa penjualan meningkat diawal bulan
2. Dengan menambakan seasionality holiday didapatkan bahwa trend peningkatan penjualan berada pada akhir tahun yaitu dibulan desember-januari

Melihat hasil perkiraan berdasarkan trend dan seonality

forecast_ts %>%
  select(ds, trend, weekly, yearly, yhat)

Model fine Tuning untuk melihat trend liner dari permintaan

Ini digunakan untuk melihat trend perkiraan dengan bentuk lain yaitu dengan menambahkan changepoint

plot(model_ts, forecast_ts) + 
  add_changepoints_to_plot(model_ts, threshold = 0) 

Garis merah putus-putus menunjukkan change point dari trend

before_2017 <- daily_demandstore03 %>% 
  mutate( 
    year = year(date) 
  ) %>% 
  filter(year < 2017) %>% 
  rename( 
    ds = "date", 
    y = "demand" 
  ) 

after_2017 <- daily_demandstore03 %>% 
  mutate( 
    year = year(date) 
  ) %>% 
  filter(year >= 2017) %>% 
  rename( 
    ds = "date", 
    y = "demand" 
  ) 
ggplot(before_2017, aes(x=ds, y=y)) + 
  geom_point() + 
  theme_minimal()

# Model sebelum 2017 dan visualisasinya
model_before_2017 <- prophet(yearly.seasonality = TRUE, 
                             changepoint.prior.scale = 0.01) %>% 
  fit.prophet(before_2017) 

## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.

future_before_2017 <- make_future_dataframe(model_before_2017, periods = 365) 
forecaset_before_2017 <- predict(model_before_2017, future_before_2017) 
plot(model_before_2017, forecaset_before_2017) + 
  add_changepoints_to_plot(model_before_2017) + 
  geom_point(data = after_2017, aes(x = as.POSIXct(ds), y=y), color = "tomato3")

# Melihat transaksi paling banyak dalam satu minggu
forecast_ts %>% 
  mutate( 
    weekday = wday(ds, label = TRUE), 
    weekly = round(weekly, 5) 
  ) %>% 
  filter(ds <= max(daily_demandstore03$date)) %>% 
  select(weekday, weekly) %>% 
  distinct() %>% 
  arrange(-weekly)

## Warning in base::check_tzones(e1, e2): 'tzone' attributes are inconsistent

# Visualisasi Model Sebelum Tahun 2017 berdasarkan hari dalam satu minggu
daily_demandstore03 %>% 
  mutate( 
    wday = wday(date, label = TRUE) 
  ) %>% 
  ggplot(aes(x=date, y=demand)) + 
  geom_point(aes(color=wday)) 

Berdasarkan visualisasi diatas dapat disumpulan bahwa penjualan tertinggi pada tiap tahun ada padari Jumat, Sabtu dan Minggu

Model Evalution

Untuk melakukan evaluasi dari model yang sudah disiapkan, maka kita bisa menguji model tersebut pada data testing.

# Menyiapkan data
test <- read.csv("data/test.csv")
glimpse(test)

## Rows: 182,500
## Columns: 4
## $ date  <chr> "2017-01-01", "2017-01-02", "2017-01-03", "2017-01-04", "2017-01~
## $ store <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1~
## $ item  <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1~
## $ sales <int> 19, 15, 10, 16, 14, 24, 14, 20, 18, 11, 14, 17, 7, 16, 29, 15, 1~

ada 182500 baris dan 4 kolom untuk data testing

# Merubah type data tanggal menjadi date
test %>% 
  mutate (date = as.Date(date)) %>%
  group_by(date) %>% 
  summarise( 
    demand=sum(sales)
  )

Bisa dilihat bahwa semua data testing ada pada tahun 2017.

Kita akan bandingkan data set Test dengan data set Training utk data sebelum tahun 2018, karena data train hanya sampai akhir 2016

cutoff <- dmy("01-01-2016") 
train <- daily_demandstore03 %>% 
  filter( 
    date < cutoff 
  ) %>% 
  rename( 
    "ds" = date, 
    "y" = demand 
  ) 
test <- daily_demandstore03 %>% 
  filter( 
    date >= cutoff 
  ) %>% 
  rename( 
    "ds" = date, 
    "y" = demand 
  ) 
ggplot(daily_demandstore03, aes(x=date, y=demand)) + 
  geom_point(data = train, aes(x=ds, y=y)) + 
  geom_point(data = test, aes(x=ds, y=y), color="tomato3")

Warna hitam adalah data training dan warna merah adalah data testing

Model Final

Model final digunakan untuk mempredisi data testing

# Menyiapkan model final
model_final <- prophet(changepoint.prior.scale = 0.05, 
                       yearly.seasonality = TRUE, 
                       holidays = holiday) %>% 
  add_seasonality(name = "monthly", period = 30.5, fourier.order = 5) %>% 
  fit.prophet(train) 

## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.

future_final <- make_future_dataframe(model_final, periods = nrow(test) + 1) 
forecast_final <- predict(model_final, future_final) 
plot(model_final, forecast_final)

Visual Model Final dengan data testing

plot(model_final, forecast_final) + 
  geom_point(data = test %>% 
               mutate(ds = as.POSIXct(ds)), aes(x=ds, y=y), color="tomato3")

Bisa dilihat bahwa model dapat memprediksi data testing dan sama dengan trend yang sudah diperkirakan sebelumnya.

Menghitung Akurasi dari Model

Untuk menentukan tingkat akurasi model, biasanya menghitung terlebih dulu nilai MAPE dari model. MAPE adalah Mean Absolute Persentase Error biasanya digunakan mengavaluasi model deret waktu.

eval <- test %>% 
  mutate( 
    ds = as.POSIXct(ds) 
  ) %>% 
  left_join(forecast_final) %>% 
  select(ds, y, yhat, yhat_upper, yhat_lower) 

## Joining, by = "ds"

eval 

Untuk mendapatkan nilai mape kita harus mengurangi nilai sebenarnya dengan nilai yang ada pada data testing (yhat) yang ada dalam data frame forecast_final

mape <- function(y, yhat) { 
  return(mean(abs(y - yhat)/ y)) 
} 
mape(eval$y, eval$yhat)

## [1] 0.04432686

Bisa dilihat bahwa nilai MAPE dari model adalah 0.04, sehingga akurasi dari model adalah 0,96 (96%). Jadi bisa disimpulkan bahwa model sangat akurat dalam melakukan prediksi.